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Abstract. We propose a direct numerical method to calculate the statistics of the 
number of transitions in stochastic processes, without having to resort to Monte Carlo 
calculations. The method is based on a generating function method, and arbitrary 
moments of the probability distribution of the number of transitions are in principle 
calculated by solving numerically a system of coupled differential equations. As an 
example, a two state model with a time-dependent transition matrix is considered and 
the first, second and third moments of the current are calculated. This calculation 
scheme is applicable for any stochastic process with a finite state space, and it would 
be helpful to study current statistics in nonequilibrium systems. 



Direct numerical method for counting statistics in stochastic processes 



2 



1. Introduction 

The statistics of nonequilibrium currents have attracted the interest of many physicists 
(for example, see [TJ.) The statistics are related to the counting of the number of specific 
transitions in stochastic processes, and these problems have been widely studied in 
physics and chemistry [2HH] . The general scheme for counting the number of transitions, 

1. e., counting statistics, has been recently developed (e.g., [5]), and various exact and 
approximate analytical results have been found. For example, a stochastic system 
under periodic perturbations can exhibit a net current, a so-called pump current, and 
it has been shown that this pump current is related to the geometric phase [9HTT]. 
For the adiabatic case in which the change of the periodic perturbations are very 
slow, an analysis with the aid of the Berry phase is valid, and analytical expressions 
for the current (the first moment) and variance have been obtained for a specific 
case [9J. However, to confirm these analytical expressions, Monte Carlo calculations 
have been performed; precisely speaking, although the first moment can be calculated 
from solutions of the master equation, Monte Carlo calculations were needed to estimate 
the variance numerically [9]. Although numerical simulations are helpful in order to 
study such time-dependent systems, Monte Carlo calculations can become very costly, 
as discussed in [9]. 

In the present paper, we derive a direct numerical method to calculate moments of 
the probability distribution of the number of specific transitions in a stochastic process. 
Although a numerical method to calculate photon emission statistics has been proposed 
recently [12], the method has only been applied to quantum systems. We will discuss 
in the following, with the aid of the generating function approach developed by Gopich 
and Szabo [5], a numerical method that applies to classical stochastic processes and 
in which the numerical effort to obtain the moments of the probability distribution 
reduces to integrating numerically a system of coupled differential equations. That is, 
this numerical method does not involve Monte Carlo type calculations. In addition, 
the method is applicable to a system under arbitrary time-dependent perturbations. It 
would be difficult in general to obtain analytic expressions for current statistics under 
complicated time-dependent perturbations, so that the direct numerical method will be 
useful for investigations of current statistics in nonequilibrium systems. 

The outline of the present paper is as follows. In section 2 we explain the general 
scheme for the generating function approach. Section 3 is the main part of the paper 
and the direct numerical method is presented by way of an example in the pump current 
problem. In addition, the validity of the proposed method is confirmed by comparison 
with the Monte Carlo simulations in section 3. Section 4 gives concluding remarks. 

2. General formalism for the generating function 

We first explain the general formalism for the generating function of the counting 
statistics. Although the formalism is similar to the one proposed in [5] , we here present 
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the general case in which the transition matrix can depend on time. 

Denoting the probability of finding the system in state n by p n (t), a master equation 
for the system is 
d 

frPnit) = E K nm (t)p m (t), (1) 

at m 

where {n n m{t)} is a transition matrix. The component (n,m) of the transition matrix 
K n m{t) is the rate of transition m — > n, and it can be time-dependent. 

We here derive the generating function for counting the number of events of a 
specific transition i A — > j\. The generalization to multiple transitions is straightforward. 

Firstly, we denote the probability, with which the system starts in state m and 
finishes in state n with N A being the number of transitions i A — > j A during time t, 
as P nm (N A \t). The probability P nm (N A \t) is obtained by repeated convolutions of the 
probability of no transitions, i.e. 

P nm (N A \t) = G' njA (t) * K jAiA (t)(? iAjA (t) • • • * K jAiA (t)(? iAjA (t) *K jAiA (t)G' iAm (t), (2) 

V v ' 

at a -i 

where g\(i) * g 2 (t) = Jq gi(t — t')g2(t')dt' denotes the convolution, and G' kl (t) is the 
probability with which the system evolves from state I to state k, provided no i A — >■ Ja 
transitions occur during time t. 

Secondly, the generating function of the probability P nm (N A \t) is defined by 

oo 

f nm (\,t)=J2^ A Pn m (N A \t). (3) 

N A 

One can see that the generating function f nm (X,t) satisfies the following integral 
equation 

f nm (X, t) = G' nm (t) + jT G' njA (t - t')\K jAiA {t')f iAm {\ t')dt'. (4) 

Thirdly, we derive the time-evolution equation for the generating function f nm (X, t). 
We notice that the probability of no transitions, G' nm (t), obeys 

^ G 'nJt) = £ KniW'Jt) - 6 n , jA K jAiA (t)G' Am (t), (5) 

where 5ij is the Kronecker delta. Hence, using the differentiation of the convolution, 

d r< , _ _ r< ( d 



dt 



jf 9l (t - t')g 2 {t')dt' = gi (0)g 2 (t) + J o l- 9l (t - f )) g 2 (t')dt', (6) 



the time-evolution equation for f nm (X,t) is derived as follows: 

9 fnm{\ t)=E Kni(t)G' im (t) ~ SnJ^j^W^t) + XG' njA (0)^ A (t)/ ixm (t) 



dt 

+ 



= J2 K ni(t)fim(^ t ) ~ $n,j A (l ~ X)Kj A i A (t)fi ATn ( X,t) , (7) 
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Figure 1. A simple stochastic model with two states. 



where we used G' nm (0) = 5 n>m . Equation (J7J) should be solved with the initial conditions 

Finally, we construct the generating function for counting the number of the target 
transitions «a — > Ja, i-e., 

oo 

F(X,t) = J2 X NA P(N A \t). (8) 

JVa=0 

Since f nm (X,t) is the generating function with specific initial state m and final state n, 
the generating function with specific final state n is constructed as 

<t> n {\t) =J2fnm{\t) Pm {0), (9) 
m 

where p m (0) is a probability distribution at initial time t = 0. Hence the generating 
function F(X,t) is calculated as 

F(A,t)=5>„(t). (10) 

n 

In addition, the time evolution equation for the generating function 4> n (X, t) is given by: 
d d 

— 4> n (X, t) = ^2 -^fnm{X, t)p m (0) 



E 



7, Knijt) fim(K t) — 8 n j A (l — X)tZj A i A (i) fi A m(X,t) 



Pm(0) 



= ^/t ni (t)0i(A,t) - 6 ndA (l - X)K jAiA (t)(j) iA (X,t), (11) 

i 

and these equations should be solved with initial conditions 0„(O) = J2 m fnm(K 0)p m (0) = 

Pn(0). 



3. Direct numerical method for the moments 



In order to explain the new direct numerical method to calculate moments in counting 
statistics, we use a simple stochastic model with two states; the discussion is easily 
extended to general cases with N states. 

The stochastic process with two states has been proposed in [H] in the context of 
the pump current problem. The system consists of three parts, as shown in Fig. [TJ 
The container can contain at most one particle. When the container is filled with one 
particle the particle can escape from the container by jumping into either one of the 
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two particle baths, i.e., the left reservoir [L] or the right one [R]. In the pump current 
problem, the transition rates k±i and k± 2 depend on time, and here we set them as 

k-i = k 2 = 1, 

ki = 1 + Rcos(ut), (12) 
fc_2 = 1 + Rsin(ut). 

Defining pi and p 2 as the respective probabilities of the container being empty or filled, 
the master equation is written as 





(13) 



/ —ki - k^ 2 k-i + k 2 
\ ki + k^ 2 —k-x - k 2 
In order to estimate the particle current from the container to the right reservoir 
[R] we need to count the number of transitions from the container to the right reservoir 
[R] and the number of transitions from [R] to the container; we should calculate the 
difference between them. According to section 2, the generating function for the particle 
current is given by 

F(A,t)=<MA,t) + 2 (A,t), (14) 
where 0i(A,t) and 02(A, t) obey the following time-evolution equations: 



d 



dt 



fe_ 2 )0i(A,*) + (k-i + M)0 2 (A,t), 
4> 2 {\,t) = (h + k-zX-^fa^t) + - k 2 )(j) 2 (\,t). 



(15) 
(16) 



Note that the transition from the container to [R] gives positive contributions to the 
statistics, so that we multiply k 2 by A. In contrast, because the transition from [R] to 
the container gives negative contributions, A -1 is multiplied to k_ 2 . 

If we obtain the explicit solution of the generating function 0j(A, t), all statistics of 
the particle current is estimated. However, in general it is difficult to obtain the explicit 
solution of the generating function. Instead of seeking analytic solutions, we proceed 
to obtain some moments directly via numerical calculations. To this effect we need to 
derive the time-evolution equations of the moments. Firstly, A in f fl5|) and ( [16]) is set to 
1: 

d 

>i|a=i + (fc-i + k 2 



dt 


di' 



>l A=l 



'2 A=l 



(-kx - k„ 2 

(h + k_ 2 )<f) 



1 1 A=l 



+ (-£;_! - k 2 



>2 A=l, 



>2 A=l) 



(17) 



where we denoted <fii(\ = 1, t) as 



A=l 



for simplicity. Because (TT71) and ( JT8l) are exactly 



the same as the original master equation (fT3|) . 



n a=i 



is interpreted as the probability 



in the original master equation, i.e., 



Pi(t) and 



>1 A=l 



p 2 (t). Next, the first 



derivatives of (|T5|) and ( !T6|) with respect to A are calculated and again A is set to 1: 



d 



dt 8X 

d d<b 2 



-ki — k_ 



90i 



dt d\ 



A=l 



A=l 



A=l 



dX 

-k- 2 p\ + (h + fc_2 



+ k 2 p 2 + (k-t + k 2 ) 



dX 



+ 



dX 
k-i - k 2 



A=l 
(902 



A=l 



dX 



(19) 
(20) 



A=l 
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Figure 2. Moments obtained by the new direct numerical method (solid lines) and 
the Monte Carlo method (filled circles with error bars), (a) first moment, (b) second 
moment, (c) third moment. 



Hence, the first moment of the particle current is calculated as follows: 



d d 



dX 



+ 



A=l 



dX 



A=lJ 



k 2 p 2 ~ fc-2Pl, 



(21) 



which recovers the well-known result [9]. 

The second and third 'factorial' moments are also estimated as follows: 



0_ 

Ot 



(n(n — 1)) 



d_ 
dt 



d 2 



OX 2 



+ 



d 2 



A=l 



dX 2 



A=l. 



2h 



dX 



+ 2k_ 2 Pi - 2/c_2 



A=l 



dX 



A=l 



(22) 



and 

d , 







_<„(„-l)(n-2)> = - 



3h 



dX 3 



+ 



d 3 c 



A=l 



dX 3 



A=l. 



dX 2 



— 6k_ 2 Pi + 6A;_ 



A=l 



dX 



3k_ 



d 2 



A=l 



dX 2 



(23) 



A=l 



The time evolution equations for d 2 <f)i/dX 2 are easily obtained, in a similar way to (|T9l) 
and ( 1201) . In numerical estimation of the time development of the differential equations, 
all d n (pi/dX n are set to initially, because 0,(A, 0) = Pi(0) does not depend on A initially. 
It is easy to obtain the second and third moments from the above factorial moments. 
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In order to calculate the m-th moment, we should solve m x 2 coupled differential 
equations. Generally, when a stochastic process consists of iV states, we need m x N 
coupled differential equations to obtain moments up to m-th order. Hence, without 
using any Monte Carlo method, we can deterministically calculate the moments for the 
particle current or the number of specific transitions. 

To confirm the validity of the new direct numerical method, we compare results 
obtained from the new method with those of the Monte Carlo method, i.e., the time- 
dependent Gillespie algorithm [13]. The parameters are R = 0.5 and u = 4.0, and the 
initial state is p\ — 1 (hence, P2 = 0); the container is empty at the initial time. In one 
data set, 10 5 Monte Carlo trajectories were used to estimate moments of the particle 
current. We repeated the Monte Carlo calculations, and collected 100 data sets. The 
error bars in Fig. [2] correspond to the standard deviation in the 100 data sets. Although 
we can calculate moments at arbitrary time in the Monte Carlo calculations, we plotted 
only some points for reference. The results of the new method and those of the Monte 
Carlo method agree completely. 

4. Concluding remarks 

In the present paper, a new numerical method to estimate moments for current statistics 
was proposed. The proposed method needs only time integrations of a system of coupled 
differential equations, so that the moments are estimated deterministically without the 
aid of Monte Carlo simulations. By way of an example, we explained the proposed 
method using a simple two-state model, but it is easy to apply this method to an 
arbitrary stochastic process with a finite state space. 

In studies of a nonequilibrium current or nonequilibrium properties, it would be 
valuable to obtain detailed information about the statistics of the current. If one 
can obtain analytical solutions for the generating function, all statistics, including the 
probability distribution for the current, are calculated from the solutions. However, 
it is difficult in general to treat a case with a complicated time-dependent transition 
matrix. On the other hand, the proposed numerical method is available to arbitrary 
time-dependency, and we believe that the proposed method is both important and useful 
in cases in which an analytical solution of the nonequilibrium current is out of reach. 
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